Discrete fourier transform in a complex vector space

ABSTRACT

An image-based phase retrieval technique has been developed that can be used on board a space based iterative transformation system. Image-based wavefront sensing is computationally demanding due to the floating-point nature of the process. The discrete Fourier transform (DFT) calculation is presented in “diagonal” form. By diagonal we mean that a transformation of basis is introduced by an application of the similarity transform of linear algebra. The current method exploits the diagonal structure of the DFT in a special way, particularly when parts of the calculation do not have to be repeated at each iteration to converge to an acceptable solution in order to focus an image.

PRIORITY CLAIM

This application claims priority under 35 U.S.C. §119 to U.S.Provisional Application Ser. No. 61/311,909 entitled “DISCRETE FOURIERTRANSFORM IN A COMPLEX VECTOR SPACE” filed on Mar. 9, 2010, the entirecontents of which are hereby incorporated by reference.

ORIGIN OF THE INVENTION

The invention described herein was made by an employee of the UnitedStates Government, and may be manufactured and used by or for theGovernment for governmental purposes without the payment of anyroyalties thereon or therefore.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This application relates to Discrete Fourier Transforms, and inparticular, to Discrete Fourier Transforms in phase retrieval for use inwaveform analysis.

2. Background

Phase retrieval is an image-based wavefront sensing method that utilizespoint-source images (or other known objects) to recover optical phaseinformation. The most famous application of phase retrieval was thediagnosis of the Hubble Space Telescope mirror edge defect discoveredsoon after the launch of Hubble and subsequent correction usingCorrective Optics Space Telescope Axial Replacement (COSTAR.) It may beironic that a phase retrieval wavefront-sensing method lies at the veryheart of the commissioning process for Hubble's successor. The earliestsuggestion of using phase retrieval as a wavefront sensing method forJWST was in 1989, nearly a year before the Hubble launch and deployment.Details documenting the Hubble Space Telescope phase retrieval analysisare discussed in the literature.

A number of image-based phase retrieval techniques have been developedthat can be classified into two general categories, the (a)iterative-transform and (b) parametric methods. Modifications to theoriginal iterative-transform approach have been based on theintroduction of a defocus diversity function or on the input-outputmethod. Various implementations of the focus-diverse iterative-transformmethod have appeared which deviate slightly by utilizing a singlewavelength or by varying the placement and number of defocused imageplanes. Modifications to the parametric approach include minimizingalternative merit functions as well as implementing a variety ofnonlinear optimization methods such as Levenburg-Marquardt, simplex, andquasi-Newton techniques. The concept behind an optical diversityfunction is to modulate a point source image in a controlled way; inprinciple, any known aberration can serve as a diversity function, butdefocus is often the simplest to implement and exhibits no angulardependence as a function of the pupil coordinates

The discrete Fourier transform (DFT) calculation is presented in“diagonal” form. By diagonal we mean that a transformation of basis isintroduced by an application of the similarity transform of linearalgebra.

What is needed is to find a more efficient implementation of the DFT forapplications using iterative transform methods, particularly phaseretrieval.

BRIEF SUMMARY

The image-based wavefront sensing method of and apparatus for focusingby determining phase difference between received images, removing thatphase difference by calculating a DFT for each signal and iterativelymultiplying the resultant until a convergence is achieved comprises thesteps of accepting point-source stellar objects to recover embeddedoptical phase information, receiving data from a sensor array whereinthe data includes embedded phase information, calculating a similaritytransform and pre-transforming the data once at the beginning of thecalculation.

The image-based wavefront sensing method of and apparatus for furtherincludes iteratively processing the data using arbitrary frequencyspacing in a diagonal fashion until a wavefront convergence is achievedto determine a phase difference by implementing the DFT in a 1dimensional linear array using just N operations, where 0<N<N×log(N)<∞,wherein N×log(N) is the number of operations of an FFT counterpart,back-transforming the data at the end of the calculation in a singlestep and focusing the received images based on the determined phase.

The general problem of diagonalizing the DFT matrix (DFT eigen-problem)has been around for years. Apparently some aspects of the problem arerelated to topics in number theory that were even solved by Gauss. Suchconsiderations have not led to more efficient implementations of the DFTover the FFT implementation, but we note that for certain applications,e.g., phase retrieval, that it is possible to exploit the diagonalstructure of the DFT in a special way, particularly when parts of thecalculation do not have to be repeated at each iteration to converge toan acceptable solution.

BRIEF DESCRIPTION OF THE SEVERAL VIEWS OF THE DRAWINGS

FIG. 1 is a functional block diagram of a system which contains adiscrete Fourier transform processor in a complex vector space,according to an example embodiment;

FIG. 2 is a flowchart representing the methodology of a discrete Fouriertransform processor in a complex vector space , according to an exampleembodiment; and

FIG. 3 is another flowchart representing the methodology of a discreteFourier transform processor in a complex vector space, according to anexample embodiment.

DETAILED DESCRIPTION

Detailed illustrative embodiments are disclosed herein. However,specific structural and functional details disclosed herein are merelyrepresentative for purposes of describing example embodiments. Exampleembodiments may, however, be embodied in many alternate forms and shouldnot be construed as limited to only the embodiments set forth herein.

Accordingly, while example embodiments are capable of variousmodifications and alternative forms, embodiments thereof are shown byway of example in the drawings and will herein be described in detail.It should be understood, however, that there is no intent to limitexample embodiments to the particular forms disclosed, but to thecontrary, example embodiments are to cover all modifications,equivalents, and alternatives falling within the scope of exampleembodiments. Like numbers refer to like elements throughout thedescription of the figures.

It will be understood that, although the terms first, second, etc. maybe used herein to describe various elements, these elements should notbe limited by these terms. These terms are only used to distinguish oneelement from another. For example, a first element could be termed asecond element, and, similarly, a second element could be termed a firstelement, without departing from the scope of example embodiments. Asused herein, the term “and/or” includes any and all combinations of oneor more of the associated listed items.

The terminology used herein is for the purpose of describing particularembodiments only and is not intended to be limiting of exampleembodiments. As used herein, the singular forms “a”, “an” and “the” areintended to include the plural forms as well, unless the context clearlyindicates otherwise. It will be further understood that the terms“comprises”, “comprising,”, “includes” and/or “including”, when usedherein, specify the presence of stated features, integers, steps,operations, elements, and/or components, but do not preclude thepresence or addition of one or more other features, integers, steps,operations, elements, components, and/or groups thereof.

Hereinafter, example embodiments will be described with reference to theattached drawings.

Example embodiments of the present invention include a Discrete FourierTransform (DFT) in phase retrieval for use in waveform analysis in aspace based application.

One such space based application under consideration to use DFT phaseretrieval methods includes the James Webb Space Telescope (JWST), one ofNASA's great observatories of the Origins program and is scheduled tolaunch later this decade. The JWST design includes using advancedoptical technology and 0-30 grade beryllium mirrors. The primary mirroris a 6.5 meter segmented primary using 18 segments in a hexagonal array.The optical design and manufacturing tolerances are specified fordiffraction-limited performance at 2 μm, and the telescope is designedfor operation over a temperature range of 30-60 K with an areal densityof less than 15 kg/m². As a result, JWST commissioning and periodicoptical maintenance incorporate san image-based wavefront sensing andcontrol (WFSC) system to align the mirror segments, minimize the effectsof figure error, and position the secondary mirror of the three-mirroranastigmat design. The wavefront sensing method specified for JWST isimage-based in the sense that point-source stellar images are collectedand processed to recover optical phase information. The primary camerafor this function is the JWST Near Infrared Camera (NIRCam), although ithas been shown that the Near Infrared Spectrograph (NIRSpec) is suitablefor this function as well.

FIG. 1 shows a block diagram of a satellite or platform in which anembodiment of the instant invention is contained. The sensor array 130is connected both physically and electronically to the platform housing120 which is in turn likewise connected the processor 110 which containsthe primary embodiments of the instant invention. The platform housingincludes all the necessary parts (not shown) of a satellite such asfuel, electrical power, communications etc. Image-based wavefrontsensing is computationally demanding due to the floating-point nature ofthe process. With certain data configurations and large array sizes, andparticularly for under-sampled data sets, the calculations can take tensof minutes or even hours to reach full convergence for the opticalwavefront. In addition, processors are generally susceptible toopto-mechanical instabilities induced by environmental factors such asjitter and turbulence. These factors not only limit fine-phasingperformance but also the stability of the phasing results. In suchcases, an obvious mitigation strategy is to increase the computationalbandwidth of the WFSC system, i.e., decrease the time required for anentire closed-loop WFSC cycle.

To improve convergence times without sacrificing accuracy and tomitigate these environmental factors, we have developed special-purposefloating-point intensive computing hardware. These processorscommunicate using a series of multi-processor computing architecturesspecifically tailored to the JWST fine-phasing methodology. Thecomputing system is based on a cluster of 64 DSPs (digital signalprocessors) operating in parallel under a shared memory architecture.The system is portable and callable from a laptop computer using atypical Ethernet connection, reaching a sustained floating-pointperformance of approximately 100 Giga-Flops with 64 DSPs fordevelopment, testing and hardware redundancy.

An embodiment of the instant invention is illustrated by some practicalaspects of calculating in a diagonal DFT basis. Aside from theapplication to iterative Fourier methods, calculating in the diagonalDFT basis also has an interesting geometrical interpretation as solvingfor the Fourier transform pair, (F_(n), f_(n)), that are co-linear or“parallel” with one another. Additional discussion on thisinterpretation is given below.

As an overview we begin by writing down the matrix form of the DFT. Nextthe hybrid-index notation of Schouten is introduced as an efficientnotation for performing the DFT calculations in a complex vector space.The similarity transformation is presented in detail for alow-dimensional DFT example (4×4 array size) and then applied toillustrate how the DFT can be carried out in a diagonal basis using thehybrid-index notation. The main results show that by “pre-transforming”the data once at the beginning of the calculation, and then“back-transforming” the data at the end of the calculation, the DFT canpotentially be implemented as a 1-d vector product using just Noperations, compared to the N×log(N) operations needed by the FFTcounterpart. The “pre-” and “back-” transforms are implemented using theDFT similarity transform and we further note that the calculation can beimplemented using arbitrary frequency spacing.

DFT as a Unitary Transform

The DFT calculation can by implemented in matrix form:

$\begin{matrix}{{{{F\left( v_{m} \right)} \equiv {{DFT}\left( v_{m} \right)}} = {\Delta \; x{\sum\limits_{n = 0}^{N - 1}{{f\left( x_{n} \right)}V_{mn}}}}},} & (1)\end{matrix}$

where V_(mn) is a non-singular M×N Vandermonde array:

V _(mn) =e ^(−i2πv) ^(m) ^((nΔx)) =e ^(−i2π(m·n)/M)=(e^(−i2π/M))^(m·n)=ω_(M) ^(m·n),   (2)

with components given by

$\begin{matrix}{V_{mn} = {\begin{pmatrix}1 & 1 & 1 & L & 1 \\1 & \omega_{M} & \omega_{M}^{2} & L & \omega_{M}^{N - 1} \\1 & \omega_{M}^{2} & \omega_{M}^{2 \cdot 2} & L & \omega_{M}^{2 \cdot {({N - 1})}} \\M & M & M & O & M \\1 & \omega_{M}^{M - 1} & \omega_{M}^{{({M - 1})} \cdot 2} & L & \omega_{M}^{{({M - 1})}{({N - 1})}}\end{pmatrix}.}} & (3)\end{matrix}$

Next we define W_(mn) as

W _(mn) =V _(mn) /√{square root over (M)}=(e ^(−i2π/M))^(m·n) /√{squareroot over (M)},   (4)

and then it is straightforward to show that W_(mn) is a unitary matrix:

W·W ^(†) =I,   (5)

using boldface matrix notation where I is the identity matrix and the“†” symbol denotes the complex-conjugate transpose or “adjoint.” Wecomment that the adjoint operation will be presented more generally inthe next Section. Next letting

f _(n) ≡Δxf(x _(n)),   (6)

Equation (1) takes the form of a unitary transformation:

$\begin{matrix}{{{{F\left( v_{m} \right)} \equiv F_{m}} = {\sum\limits_{n = 0}^{N - 1}{W_{mn}f_{n}}}},} & (7)\end{matrix}$

with the inverse transformation given by (generally but not always wetake M=N):

$\begin{matrix}{f_{n} = {\sum\limits_{m = 0}^{M - 1}{\left( W_{mn} \right)^{\dagger}{F_{m}.}}}} & (8)\end{matrix}$

Given (5) and (7), the DFT can be analyzed using the transformationtheory of complex vector spaces. The hybrid-index notation of Schoutensupplies a natural framework for this approach and is presented in thenext Section.

Hybrid Index Notation

As illustrated above, the DFT is a unitary transformation in a complexvector space. Therefore, the calculations can be performed elegantly andefficiently using the tensor-transformation formalism of Schouten. Inparticular, Schouten's hybrid-index notation provides a naturalframework for the DFT calculations. The goal in using this approach isto give a more general representation of the calculations by identifyingcomponents with quantities of a geometric manifold. Often, but notalways, such identification can lead to alternative or more efficientcalculation strategies as discussed later. Therefore, the goal in thenext Section is to sort out the details of the hybrid-index notation forthe DFT, and thus identify the DFT with quantities in a geometricmanifold. In addition, by adopting a notation that makes thetransformation rules of a quantity explicit, the meaning of thecalculations can often be clarified and the notation serves as anefficient bookkeeping device for checking consistency throughout thecalculations.

DFT in Hybrid-Index Notation

The DFT is defined here as a unitary transform between conjugate Fourierdomains. The 1-d DFT transform is then expressed using Schouten'shybrid-index notation as

F ^(m) =W ^(m) _(·n) f ^(n),   (9)

using the Einstein sum convention defined over repeated indices in anupper/lower combination:

$\begin{matrix}{{\sum\limits_{n}^{\;}{W_{mn}f_{n}}} \equiv {W_{\cdot n}^{m}{f^{n}.}}} & (10)\end{matrix}$

The lower/upper indices denote the covariant and contravarianttransformation rules—the upper indices are not exponents in thisnotation. Exponents of ( a quantity will be denoted using parentheses,or brackets, respectively, for either a transformation of coordinates,or a “point transformation.” Allowable transformations of coordinatesfor this application will be discussed later. The distinction betweenthe transformations of coordinates and point transformations arefundamental to Schouten's kernel-index method. To elaborate briefly,coordinate transformations are denoted by a transformation of basis,say, m→m′, where the new basis set is denoted by m′. The components inthe new basis are thus

u ^(m) ′=A _(m) ^(m) ′u _(m),   (11)

where A_(m) ^(M)′ is the transformation matrix and the u^(m)′ are thevector components in the new basis set (u is arbitrary).

A point transformation, on the other hand, is denoted by

U ^(m) =P _(·n) ^(M) u ^(n),   (12)

and is the form adopted for the DFT in (9). Another importantdistinction between a transformations of basis and a pointtransformation is that the indices for the point transformation areordered and the positions are marked using over- or under-dots as inP_(·n) ^(m) (12). The main take-away is that for coordinatetransformations, the kernel symbol, u, remains the same as in (11). Butfor point transformations, the kernel symbol changes (as in (12)) butthe indices remain the same. So initially at least in this context, theDFT transformation in (9) is defined as a point transformation, butlater we find that when calculating the DFT in a diagonal basis, it isnatural to reconsider (9) using a common kernel symbol for the Fouriertransform pair, specifically, (F^(m), f^(n))→(f _(m) , f^(n)). We notealso that it is possible to mix the two transformations (basis andpoint) as is required when considering the allowable transformations ofbasis in digital signal processing. These details will be discussed inmore detail later.

Continuing with an explanation of the formalism, the complex conjugateof (9) is given by

F ^(m) = W ^(m) _(n) f ^(n) ≡ f ^(m,)  (13)

using the over-bar to denote the complex conjugation of components andkernel symbols in accordance with Schouten's hybrid-index notation andnoting for consistency that

F ^(m) = F ^(m)

F^(m) =

=F ^(m).   (14)

The conjugation of indices and kernel objects is fundamental to thestudy of transformations in a complex vector space, and as illustratedbelow, is perfectly suited for calculations of the DFT transform.

We choose an initial convention that the data values are labeled bycontravariant components, f^(m), and then by complex conjugation of thekernel symbol and indices in various combinations, several other typesof components are represented in Schouten's notation by the “raising”and “lowering” of indices using the “fundamental” or metric tensor, α_(kl):

f^(m),f _(m) , f ^(m) , f _(m),   (15)

where the latter two follow by complex conjugation of the first two.Specifically, we mean that

f ^(m)

f _(m) =a _(mm) f ^(m) ; f ^(m)

f _(m) =a _(mm) f ^(m) ,   (16)

As mentioned above, co- and contravariant transformation rules can beswitched back and forth using the metric or fundamental tensor, a _(kl),to raise or lower indices as appropriate. The metric tensor also defineshow the differential element of distance is calculated in the manifoldas well as the inner product between two vectors, say, u^(k) and v^(l):

u|v

=a _(kl) ū ^(k) v ^(l) =ū _(l) v ^(t) =ū ^(k) v _(k) ,   (17)

pointing out that in a linear complex vector space, the “bra-ket”notation of Dirac is easily incorporated in this approach as illustratedon the LHS of (17). Illustrating the use of the metric for raising andlowering indices we have:

W _(·n) ^(m) =a _(ln) W ^(m l) and f ^(n) =ā ^(n k) f _(k) ,   (18)

and then substituting these expressions back into (9) we see that analternative form of the DFT is obtained:

F ^(m) =W _(·n) ^(m) f ^(n)=(a _(ln) W ^(m l) )(a ^(n k) f _(k) )=a_(ln) a ^(n k) W ^(m l) f _(k) =δ _(l) ^(· k) W ^(m l) f _(k)

F ^(m) =W ^(m n) f _(n) .   (19)

Adjoint in a Geometric Manifold

In general, the adjoint operation applied to the W_(·n) ^(m), denoted by{tilde over (W)}_(·n) ^(m)=(W_(·n) ^(m))^(†), is defined through atwo-step operation using both the metric tensor and complex conjugation:

-   -   (a) each index is raised or lowered as appropriate using the        metric:

a _(nk) a ^(l m) W _(·l) ^(k) =W _(n) ^(· m) ,   (20)

-   -   (b) the result in (a) is then complex conjugated to give

(W _(·n) ^(m))^(†=)

⁼

⁼ W _(n) ^(·m), tm (21)

and finally we write

{tilde over (W)} _(·n) ^(m)=≡(W _(·n) ^(m))^(†) = W _(n) ^(·m),   (22)

For comparison, in matrix boldface notation:

(W)^(\) = W ^(T) =a Wa ⁻¹ =W ³¹ ¹,   (23)

where “T” defines the usual transpose operation. We note especially thatthe example in (23) illustrates a degree of non-generality in the waythe adjoint is usually calculated in matrix notation. Specifically,defining the adjoint using components in a geometric manifold (usingsteps (a) and (b) above) can give different results than simply complexconjugating and then performing the transpose operation using the “T” asin (23). In special coordinates, a _(nk)≐δ _(nk), these two areidentical, but in general, that is not always the case for arbitrary a_(nk).

Since the W_(·n) ^(m) are unitary, the inverse,

${\overset{- 1}{W}}_{\cdot n}^{m},$

is defined using the adjoint and is expressed (from (22)):

$\begin{matrix}{{{\overset{\sim}{W}}_{\cdot n}^{m} = {{\overset{- 1}{W}}_{\cdot n}^{m} = {\left( W_{\cdot n}^{m} \right)^{\dagger} = {\overset{\_}{W}}_{n}^{\cdot m}}}},} & (24)\end{matrix}$

such that

W _(·k) ^(n)

=δ_(·k) ^(l).  (25)

We further comment that Schouten did not initially use an over-bar onthe kernel symbol as in we have in (25). Considering the W_(·n) ^(m)using two upper indices, W^(m n) (see (18) and (19)) the correspondingadjoint is similarly calculated:

(W ^(m n) )^(\) =a _(mk) a _(ln) W ^(k l) =

= W _(m n) =

.  (26)

The inverse DFT transformation to (9) is thus defined:

{tilde over (W)} _(·m) ^(n) F ^(m)=

W _(·k) ^(m) f ^(k)δ_(k) ^(n) f ^(k) =f ^(n),   (27)

as expected, or when using the form of (19) the inverse DFT is expressed

{tilde over (W)} _(m n) F ^(m)=

(W ^(m k) f _(k) )=δ _(n) ^(k) f _(k) =f _(n) .   (28)

A few additional comments from this Section: in tensor transformationtheory, there is geometric meaning attributed to these differing butcomplementary transformation rules. Specific transformations are definedby an allowable Group of transformations, G_(a), as discussed further inSchouten. The corresponding allowable transformations of basis for theDFT application (digital signal processing) can include fractional andinteger shifts of the data origin in addition to scaling and phasetransformations. The details for these various transformations in thisnotation will be discussed in a later set of notes.

Matrix Multiplication Using Index Notation

When calculating with components, index order is important for propermatrix multiplication and thus we adopt the following convention for arank-2 quantity:

$\begin{matrix}{A_{\cdot \underset{\underset{col}{\downarrow}}{n}}^{m\rightarrow{row}},} & (29)\end{matrix}$

and the inverse of A_(·n) ^(m) is notated as,

${\overset{- 1}{A}}_{\cdot n}^{m},$

such that

${A_{\cdot n}^{m}{\overset{- 1}{A}}_{\cdot k}^{n}} = \delta_{\cdot k}^{m}$

where δ_(·k) ^(m) is the Kronecker delta. To illustrate, the matrixmultiplication of two matrices A and B is notated as:

C _(·k) ^(m) =A _(·n) ^(m) B _(·k) ^(n), (30)

and to be specific we calculate (30) using Mathematica with 2×2 arrays:

Table [Sum [(A[[um,n]]B[[n,lk]]), {n,1,NN}], {um,1,2}, {lk,1, 2},]]  (31)

where um denotes “upper-m” and lk denotes “lower-k;” the upper/lowerdummy variable, n, is left un-notated in Mathematica with regard to itsposition. The result is of course (30):

$\begin{matrix}{{C_{\cdot k}^{m} = {{A_{\cdot n}^{m}B_{\cdot k}^{n}} = {{AB} = \begin{pmatrix}{{a_{\cdot 1}^{1}b_{\cdot 1}^{1}} + {a_{\cdot 2}^{1}b_{\cdot 1}^{2}}} & {{a_{\cdot 1}^{1}b_{\cdot 2}^{1}} + {a_{\cdot 2}^{1}b_{\cdot 2}^{2}}} \\{{a_{\cdot 1}^{2}b_{\cdot 1}^{1}} + {a_{\cdot 2}^{2}b_{\cdot 1}^{2}}} & {{a_{\cdot 1}^{2}b_{\cdot 2}^{1}} + {a_{\cdot 2}^{2}b_{\cdot 2}^{2}}}\end{pmatrix}}}},} & (32)\end{matrix}$

where AB is the conventional matrix product using matrix boldfacenotation. For comparison, transposing both the m-n and n-k indexpositions in (30) reverses the order and gives the BA matrix product:

$\begin{matrix}{{C_{\cdot k}^{m} = {{A_{n}^{\cdot m}B_{k}^{\cdot n}} = {{BA} = \begin{pmatrix}{{b_{\cdot 1}^{1}a_{\cdot 1}^{1}} + {b_{\cdot 2}^{1}a_{\cdot 1}^{2}}} & {{b_{\cdot 1}^{1}a_{\cdot 2}^{1}} + {b_{\cdot 2}^{1}a_{\cdot 2}^{2}}} \\{{b_{\cdot 1}^{2}a_{\cdot 1}^{1}} + {b_{\cdot 2}^{2}a_{\cdot 1}^{2}}} & {{b_{\cdot 1}^{2}a_{\cdot 2}^{1}} + {b_{\cdot 2}^{2}a_{\cdot 2}^{2}}}\end{pmatrix}}}},} & (33)\end{matrix}$

and the corresponding Mathematica code for (33) is

Table[Sum[(A[[n,um]]B[[lk,n]]), {n,1,NN}], {um,1,2}, {lk,1,2}]  (34)

Metric Tensor Deriving from Unitary Invariance

In summary of Schouten's notation, various geometric quantities aredefined that are distinguished by their transformation rules using indexposition. A group of allowable transformations, G_(a), is introduced andthen several invariant quantities may be constructed based on invarianceof a given symmetry group. One such quantity is the metric tensor, a_(kl), defining “distance” in the manifold and also defining the innerproduct between any two vectors. Therefore, one goal in these notes isto solve for the DFT metric tensor, a _(kl), by postulating nothing morethan unitary invariance under the action of the point transformation,W_(·n) ^(m), as illustrated below. Specifically, given that the W_(·n)^(m) are defined as a unitary transformation, we introduce a metrictensor defined by invariance under the action of the W_(·n) ^(m)transformation such that

W _(·l) ^(m) W _(k) ^(· n) a _(nm) =a _(kl),   (35)

or in matrix form:

W a W ⁻¹ =a.  (36)

Quantities that are invariant under the action of the unitary group arehermitian and thus the a _(kl) satisfy

a _(kl)=(a kl)^(\) =ā _(lk),   (37)

The introduction of a hermitian metric a _(kl) into the E_(n) geometricmanifold defines the manifold as a U_(n) which differs from the earliernotation used by Schouten where this manifold was labeled as the {tildeover (R)}_(n).

General form of Parseval's Theorem

We can show that Parseval's theorem is really implied by unitaryinvariance, and not necessarily “energy conservation” as it is commonlyreferred to in the engineering literature. To begin, we consider themagnitude squared of F^(m)→|F|², giving

F|F

=|F| ² =a _(mn) F ^(m) F ^(n).  (38)

Substituting the DFT transformation rule from (9) into (38) we have

|F| ² =a _(nm) W _(·k) ^(m) W _(l) ^(· n) f ^(k) f ^(l) ,   (39)

but since the a _(k) _(l) are invariant under this unitarytransformation (see (35)) then the RHS of (39) simplifies and we get

|F| ² =a _(kl) f ^(k) f ^(l).   (40)

Therefore, by comparison with (38), Parseval's theorem is seen to holdfor any solution of the a _(k) _(l)as implied by unitary invariance(35), and thus in general:

a _(mn) F ^(n) F ^(n) =a _(kl) f ^(k) f ^(l).   (41)

The common form of Parseval's theorem is expressed in vector notation as

|F| ² =|f| ².  (42)

To give some perspective on these results we look forward a bit andcomment that the standard DFT approach adopts a special solution(coordinate system) for the a _(kl) given by the Kronecker delta (asshown later):

a _(nk)≐δ _(nk)  (43)

and so in index notation (42) takes the familiar “Euclidean” form:

δ _(mn) F ^(m) F ^(n)=δ _(kl) f ^(k) f ^(l),   (44)

or as it is commonly written:

|F ₁|² +|F ₂|² +L+|F _(N)|² =|f ₁|² +|f ₂|² +L+|f _(N)|².  (45)

The point is that for a general a _(kl), Parseval's theorem (41) willnot necessarily take the Euclidean form as demonstrated later below (see(67)), although unitary invariance implies that (41) still holds.

A couple of additional comments: as discussed later the δ _(kl) is asolution for the a _(kl) and corresponds to a positive definitesignature solution. The approach taken here is motivated by that factthat in general, the eigen-problem requires the use of a general metric,and by no means does this necessarily correspond to the Kronecker deltain a unitary complex vector space, U_(n). But even so, we later adoptthe Kronecker delta solution to make progress on the main point ofcalculating the DFT in a diagonal basis, and for easy comparison withpreviously published results.

Metric Solutions

By postulating unitary invariance of the metric we can solve for the a_(kl) using (35) and (37). The W_(·k) ^(m) are known throughidentification with (4):

$\begin{matrix}{{W_{\cdot k}^{m} \equiv {\frac{1}{\sqrt{M}}\left( ^{{- {{2\pi}{({mk})}}}/M} \right)}} = {\frac{1}{\sqrt{2}}{\begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}.}}} & (46)\end{matrix}$

As discussed above, the inverse transformation of (46) is definedthrough the adjoint operation in (20). But at this point in thediscussion we have not yet solved for the metric a _(kl), which isneeded to form the adjoint. Therefore, we first solve for the

${\overset{- 1}{W}}_{\cdot k}^{n}$

such that

$\begin{matrix}{{{W_{\cdot n}^{m}{\overset{- 1}{W}}_{\cdot k}^{n}} = \delta_{\cdot k}^{m}},} & (47)\end{matrix}$

and the

${\overset{- 1}{W}}_{\cdot k}^{n}$

are given by:

$\begin{matrix}{{\overset{- 1}{W}}_{\cdot k}^{m} = {{\overset{\_}{W}}_{k}^{\cdot m} = {{\frac{1}{\sqrt{M}}\left( ^{{{2\pi}{({mk})}}/M} \right)} = {\frac{1}{\sqrt{2}}{\begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}.}}}}} & (48)\end{matrix}$

Continuing to the general

×2 metric, a _(kl)

defined as

$\begin{matrix}{{a_{\overset{\_}{k}l} = {\begin{pmatrix}z_{\overset{\_}{1}1} & z_{\overset{\_}{1}2} \\z_{\overset{\_}{2}1} & z_{\overset{\_}{2}2}\end{pmatrix} = \begin{pmatrix}{x_{11} + {\; y_{11}}} & {x_{12} + {\; y_{12}}} \\{x_{21} + {\; y_{21}}} & {x_{22} + {\; y_{22}}}\end{pmatrix}}},} & (49)\end{matrix}$

where the complex notation is dropped on indices of real components. Thecondition for unitary invariance (35) then gives the followingsimultaneous system of 4 equations and 8 unknowns:

$\begin{matrix}{{\frac{1}{2}\begin{pmatrix}{z_{\overset{\_}{1}1} + z_{\overset{\_}{1}2} + z_{\overset{\_}{2}1} + z_{\overset{\_}{2}2}} & {z_{\overset{\_}{1}1} - z_{\overset{\_}{1}2} + z_{\overset{\_}{2}1} - z_{\overset{\_}{2}2}} \\{z_{\overset{\_}{1}1} + z_{\overset{\_}{1}2} - z_{\overset{\_}{2}1} - z_{\overset{\_}{2}2}} & {z_{\overset{\_}{1}1} - z_{\overset{\_}{1}2} - z_{\overset{\_}{2}1} + z_{\overset{\_}{2}2}}\end{pmatrix}} = {\begin{pmatrix}z_{\overset{\_}{1}1} & z_{\overset{\_}{1}2} \\z_{\overset{\_}{2}1} & z_{\overset{\_}{2}2}\end{pmatrix}.}} & (50)\end{matrix}$

The system (50) is under-determined and can therefore have many possiblesolutions but will be further constrained by the hermitian condition(37). Solving (50) eliminates 4 of the unknowns (z ₁₁ and z ₁₂) in termsof the other 4 unknowns (z ₂₁ and z ₂₂) showing that a _(kl) has thegeneral form:

$\begin{matrix}{a_{\overset{\_}{k}l} = {\begin{pmatrix}{{2z_{\overset{\_}{2}1}} + z_{\overset{\_}{2}2}} & z_{\overset{\_}{2}1} \\z_{\overset{\_}{2}1} & z_{\overset{\_}{2}2}\end{pmatrix}.}} & (51)\end{matrix}$

Next applying the hermitian (symmetry) condition (35) to (51) gives thecondition that the 2×2 a _(kl) must be real valued:

$\begin{matrix}{{\begin{pmatrix}{{2y_{21}} + y_{22}} & y_{21} \\y_{21} & y_{22}\end{pmatrix} = {\left. \begin{pmatrix}0 & 0 \\0 & 0\end{pmatrix}\Rightarrow{{Im}\left( a_{\overset{\_}{k}l} \right)} \right. = 0}},} & (52)\end{matrix}$

but as expected there are multiple solutions for the a _(kl) in terms ofthe 2 final parameters x₂₁ and x₂₂. The general solution for the 2×2case is thus reduced from 8 to 2 unspecified real components and isgiven by

$\begin{matrix}{a_{\overset{\_}{k}l} = {\begin{pmatrix}{{2x_{21}} + x_{22}} & x_{21} \\x_{21} & x_{22}\end{pmatrix}.}} & (53)\end{matrix}$

As a check, substituting (53) back into (35) verifies that this a _(kl)is invariant under the action of the W_(·k) ^(m) for arbitrary (butreal) numerical values of x₂₁ and x₂₂ in (53).

We have shown that an infinite number of solutions exist for the 2×2metric, a _(kl), as implied by unitary invariance in terms of two realparameters. We note that one of the simplest and most importantsolutions is a special case of (53) when x₂₁=0 and x₂₂=1 giving theKronecker delta as mentioned earlier:

$\begin{matrix}{{{{\underset{1}{a}}_{\overset{\_}{k}l}_{\underset{x_{22}\rightarrow 1}{x_{21}\rightarrow 0}}}\overset{*}{=}{\begin{pmatrix}1 & 0 \\0 & 1\end{pmatrix}\overset{*}{=}\delta_{{\overset{\_}{k}l}\;}}},} & (54)\end{matrix}$

where the complex conjugation is kept for consistency, and the numeral“1” beneath the kernel symbol in Schouten's notation labels this case asa distinct solution, different from other solutions, say

${\underset{2}{a}}_{\overset{\_}{k}l}.$

The labeling of kernel symbols in this way is used extensively in theeigen-problem discussed later.

We comment that since (54) is positive definite, as discussed at greaterlength below in the Section titled Generalized DFT, there will be nodistinction between the other various quantities that may be derivedfrom the a _(kl) in the U₂, and thus, no numerical distinction existsfor components labeled by upper/lower indices, or even the complexconjugation of indices. This is seen in the solution (54) above.

As illustrated above there are an infinite number of solutions for the a_(kl). Several specific solutions are shown in (55) based on variouschoices for the x₂₁ and x₂₂ components:

$\begin{matrix}{{{{{\underset{1}{a}}_{\overset{\_}{k}l}_{\underset{x_{22}\rightarrow 1}{x_{21}\rightarrow 0}}}\overset{*}{=}\begin{pmatrix}1 & 0 \\0 & 1\end{pmatrix}};{{{\underset{2}{a}}_{\overset{\_}{k}l}_{\underset{x_{22}\rightarrow{- 1}}{x_{21}\rightarrow 1}}}\overset{*}{=}\begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}}}{{{{{\underset{3}{a}}_{\overset{\_}{k}l}_{\underset{x_{22}\rightarrow 0}{x_{21}\rightarrow 1}}}\overset{*}{=}\begin{pmatrix}2 & 1 \\1 & 0\end{pmatrix}};{{{\underset{4}{a}}_{\overset{\_}{k}l}_{\underset{x_{22}\rightarrow 1}{x_{21}\rightarrow{- \frac{1}{2}}}}}\overset{*}{=}\begin{pmatrix}0 & {- \frac{1}{2}} \\{- \frac{1}{2}} & 1\end{pmatrix}}},}} & (55)\end{matrix}$

Given some specific solutions for the a _(kl) in (55), we examinevarious other quantities that are derived in Schouten's notation. Forinstance, consider the second solution in (55) that is a multiple of theW_(·n) ^(m):

$\begin{matrix}{{{a_{\overset{\_}{k}l} \equiv {\underset{2}{a}}_{\overset{\_}{k}l}}\overset{*}{=}{\left. \begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}\Rightarrow{\overset{\_}{a}}^{k\overset{\_}{l}} \right.\overset{*}{=}{\frac{1}{2}\begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}}}},} & (56)\end{matrix}$

and then for generality define a Complex data vector:

f ^(n)=(z ¹ ,z ²)=(x ¹ +iy ¹ ,x ² +iy ²)  (57)

The DFT transformation (9) is thus given by

$\begin{matrix}{F^{m} = {{W_{\cdot n}^{m}f^{n}} = {\frac{1}{\sqrt{2}}{\left( {{z^{1} + z^{2}},{z^{1} - z^{2}}} \right).}}}} & (58)\end{matrix}$

As a consistency check on both the solution for the a _(kl) of (56), andthe legitimacy of constructing the inverse transformation by use of theadjoint operation in (20), i.e.,

$\mspace{79mu} {{{\overset{- 1}{W}}_{\cdot n}^{m} = {\left( W_{\cdot n}^{m} \right)^{\dagger} \equiv \text{?}}},{\text{?}\text{indicates text missing or illegible when filed}}}$

we can check the inverse transformation calculated by way of theadjoint:

$\begin{matrix}{\mspace{79mu} {\begin{matrix}{\left( W_{\cdot n}^{m} \right)^{\dagger} = {a_{\overset{\_}{n}k}a^{l\overset{\_}{m}}W_{\cdot l}^{k}}} \\{= {\text{?}\frac{1}{\sqrt{2}}\text{?}\frac{1}{2}\text{?}}} \\{= {\frac{1}{\sqrt{2}}\begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}}} \\{{= {\overset{- 1}{W}}_{\cdot n}^{m}},}\end{matrix}{\text{?}\text{indicates text missing or illegible when filed}}}} & (59)\end{matrix}$

in agreement with (48). The inverse

transformation is also

for consistency giving

$\begin{matrix}{{{\overset{\sim}{W}}_{\cdot n}^{m}F^{m}} = {{\frac{1}{\sqrt{2}}\begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}\begin{pmatrix}{z^{1} + z^{2}} \\{z^{1} - z^{2}}\end{pmatrix}} = {\left( {z^{1},z^{2}} \right) = {f^{n}.}}}} & (60)\end{matrix}$

Several other quantity are derived from the f^(n) using the a _(kl):

f ^(n)=(z ¹ ,z ²); f _(n) =a _(nk) f ^(k)=(z ¹ +z ² ,z ¹ −z ²)

= f ^(k) =( z ¹ , z ² ); f _(k) =a _(nk) f ^(n) =( z ¹ + z ² , z ¹ − z ²),  (61)

and comparing the RHS's in (61), as expected we have f _(n) =

, showing consistency using the hybrid-index notation. Correspondingly,various other F^(m) quantities are given by

$\begin{matrix}{{{F^{m} = {\frac{1}{\sqrt{2}}\left( {{z^{1} + z^{2}},{z^{1} - z^{2}}} \right)}};{F_{\overset{\_}{m}} = {{a_{\overset{\_}{m}k}F^{k}} = {\sqrt{2}\left( {z^{1},z^{2}} \right)}}}}{{{\overset{\_}{F}}^{\overset{\_}{m}} = {\frac{1}{\sqrt{2}}\left( {{{\overset{\_}{z}}^{\overset{\_}{1}} + {\overset{\_}{z}}^{\overset{\_}{2}}},{{\overset{\_}{z}}^{\overset{\_}{1}} - {\overset{\_}{z}}^{\overset{\_}{2}}}} \right)}};{{\overset{\_}{F}}_{m} = {{a_{\overset{\_}{n}m}{\overset{\_}{F}}^{\overset{\_}{n}}} = {\sqrt{2}\left( {{\overset{\_}{z}}^{\overset{\_}{1}},{\overset{\_}{z}}^{\overset{\_}{2}}} \right)}}}}} & (62)\end{matrix}$

and components for the W^(m n) and the W _(mn) are also calculated:

$\begin{matrix}{{{W^{m\overset{\_}{n}} = {{a^{k\overset{\_}{n}}W_{\cdot k}^{m}} = \begin{pmatrix}\frac{1}{\sqrt{2}} & 0 \\0 & \frac{1}{\sqrt{2}}\end{pmatrix}}};}{W_{\overset{\_}{m}n} = {{a_{\overset{\_}{k}n}W_{\overset{\_}{m}}^{\cdot \overset{\_}{k}}} = {\begin{pmatrix}\frac{1}{\sqrt{2}} & 0 \\0 & \frac{1}{\sqrt{2}}\end{pmatrix}.}}}} & (63)\end{matrix}$

Using the f _(n) from (61) and then comparing with the F^(m) in (62) wecan check the forward DFT for consistency using the W^(m n) of (63):

$\begin{matrix}\begin{matrix}{{W^{m\overset{\_}{n}}f_{\overset{\_}{n}}} = {\begin{pmatrix}\frac{1}{\sqrt{2}} & 0 \\0 & \frac{1}{\sqrt{2}}\end{pmatrix}\begin{pmatrix}{z^{1} + z^{2}} \\{z^{1} - z^{2}}\end{pmatrix}}} \\{= {\frac{1}{\sqrt{2}}\left( {{z^{1} + z^{2}},{z^{1} - z^{2}}} \right)}} \\{= {F^{m}.}}\end{matrix} & (64)\end{matrix}$

The adjoint operation on the W^(m n) and W _(mn) of (63) is similarlyverified:

$\begin{matrix}{\mspace{79mu} {\begin{matrix}{\left( W_{\cdot n}^{m} \right)^{\dagger} = \text{?}} \\{= {a_{m\overset{\_}{k}}a_{\overset{\_}{n}l}W^{m\overset{\_}{n}}}} \\{= {\begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}\begin{pmatrix}\frac{1}{\sqrt{2}} & 0 \\0 & \frac{1}{\sqrt{2}}\end{pmatrix}\frac{1}{2}\begin{pmatrix}1 & 1 \\1 & {- 1}\end{pmatrix}}} \\{{= \begin{pmatrix}\sqrt{2} & 0 \\0 & \sqrt{2}\end{pmatrix}},}\end{matrix}{\text{?}\text{indicates text missing or illegible when filed}}}} & (65)\end{matrix}$

such that

W ^(m n)

=δ _(k) ^(n) .   (66)

Finally, Parseval's theorem is demonstrated:

a _(mn) F ^(m) F ^(n) =a _(kl) f ^(k) f ^(l) =|z ¹|² −|z ²|^(s),   (67)

which is not in Euclidean form but is nevertheless invariant.

Comments on a General DFT

Further restrictions may be placed on the solution values for the x₂₁and x₂₂ in (53), in accordance with the classification of the a _(kl) asbeing positive definite, negative definite, or indefinite, correspondingto the number of values (or signature)>0, <0, or mixed, respectively,along the diagonal of the a _(kl)). Note that when a _(kl) is definitethat there is no distinction between the various quantities that may bederived from a _(kl) in the U₂ and thus, no numerical distinction existsfor components labeled by upper/lower indices, or even the complexconjugation of indices; these notations are merely formal in thisinstance. We comment further that since the a _(kl) are required whensolving the eigen-problem in general (more details on this below) therecan in principle at least, be different solutions for the eigen-problemcorresponding to different solutions for the a _(kl). Note that allearlier work on the DFT eigen-problem uses the a _(kl)≐δ _(kl) solutionwithout considering the possibility that the a _(kl) might be defined ina more general way based on unitary invariance. Thus, differentsolutions for the a _(kl) may lead to possibly different DFTeigen-quantities, but all the while preserving unitary invariance andthe hermitian character of the a _(kl). Thus, from the perspective ofthe transformation theory of complex vector spaces, we can regard someproperties of the DFT as deriving from the characteristics of the a_(kl), i.e., that (a) a _(kl) is invariant under unitary transformationsand (b) is hermitian. Condition (b) really derives from (a) but it isconvenient when solving for the a _(kl) components (as illustratedabove) to regard (b) as an independent condition.

Similarly as in the 2×2 example, we can solve for the a _(kl)corresponding to higher dimensional data sets. Of special interest aresolutions in which the a _(kl) are complex valued since these solutionsprovide further consistency checks on the proper use of Schouten's'notation. We proceed by deriving the 4×4 W_(·k) ^(m) and its inverse,

${\overset{- 1}{W}}_{\cdot k}^{m},$

from (46):

$\begin{matrix}{{{W_{\cdot k}^{m} = {\frac{1}{2}\begin{pmatrix}1 & 1 & 1 & 1 \\1 & {- i} & {- 1} & i \\1 & {- 1} & 1 & {- 1} \\1 & i & {- 1} & {- i}\end{pmatrix}}};}\begin{matrix}{{\overset{- 1}{W}}_{\cdot k}^{m} = \text{?}} \\{= {\overset{\_}{W}}_{k}^{\cdot m}} \\{{= {\frac{1}{2}\begin{pmatrix}1 & 1 & 1 & 1 \\1 & i & {- 1} & {- i} \\1 & {- 1} & 1 & {- 1} \\1 & {- i} & {- 1} & i\end{pmatrix}}},}\end{matrix}{\text{?}\text{indicates text missing or illegible when filed}}} & (68)\end{matrix}$

and a general 4×4 metric is defined by

$\begin{matrix}{a_{\overset{\_}{k}l}\begin{pmatrix}z_{\overset{\_}{1}1} & z_{\overset{\_}{1}2} & z_{\overset{\_}{1}3} & z_{\overset{\_}{1}4} \\z_{\overset{\_}{2}1} & z_{\overset{\_}{2}2} & z_{\overset{\_}{2}3} & z_{\overset{\_}{2}4} \\z_{\overset{\_}{3}1} & z_{\overset{\_}{3}2} & z_{\overset{\_}{3}3} & z_{\overset{\_}{3}4} \\z_{\overset{\_}{4}1} & z_{\overset{\_}{4}2} & z_{\overset{\_}{4}3} & z_{\overset{\_}{4}4}\end{pmatrix}} & (69)\end{matrix}$

The condition for unitary invariance (35) then gives a simultaneoussystem of 16 equations and 32 unknowns which can be further reduced to 6free parameters by imposing the hermitian condition (37) on the a _(kl).After considerable algebra the general solution to (35) for the a _(kl),constrained by the hermitian condition (37) is given by

$\begin{matrix}{a_{\overset{\_}{k}l} = {\begin{pmatrix}{x_{33} + {2\left( {x_{41} + x_{43}} \right)}} & {x_{41} - {iy}_{43}} & \left( {x_{33} + x_{41} - x_{42} + x_{43} - x_{44} + {2{iy}_{43}}} \right) & {x_{41} - {iy}_{43}} \\{x_{41} + {iy}_{43}} & x_{44} & {x_{43} + {iy}_{43}} & x_{42} \\\left( {x_{33} + x_{41} - x_{42} + x_{43} - x_{44} + {2{iy}_{43}}} \right) & {x_{43} - {iy}_{43}} & x_{33} & {x_{43} - {iy}_{43}} \\{x_{43} + {iy}_{43}} & x_{42} & {x_{43} + {iy}_{43}} & x_{44}\end{pmatrix}.}} & (70)\end{matrix}$

As before, there are an infinite number of both real and complex valued4×4 solutions for the a _(kl) in terms of the 6 free parameters in (70).Several example solutions corresponding to basic numerical values forthe x and y components are listed in (71):

$\begin{matrix}{{\underset{1}{a}}_{k\overset{\_}{l}}{_{\underset{{x_{43}\rightarrow 0},{y_{43}\rightarrow 0}}{\underset{{x_{41}\rightarrow 0},{x_{42}\rightarrow 0}}{{x_{33}\rightarrow 1},{x_{44}\rightarrow 1}}}}{{\overset{*}{=}\begin{pmatrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{pmatrix}};{{\underset{2}{a}}_{k\overset{\_}{l}}{_{\underset{{x_{43}\rightarrow 0},{y_{43}\rightarrow 1}}{\underset{{x_{41}\rightarrow 0},{x_{42}\rightarrow 0}}{{x_{33}\rightarrow 1},{x_{44}\rightarrow 1}}}}{\overset{*}{=}{\begin{pmatrix}1 & {- i} & {2i} & {- i} \\i & 1 & i & 0 \\{{- 2}i} & {- i} & 1 & {- i} \\i & 0 & i & 1\end{pmatrix}{\underset{3}{a}}_{k\overset{\_}{l}}{_{\underset{{x_{43}\rightarrow 1},{y_{43}\rightarrow 0}}{\underset{x_{41,}\rightarrow{{{- 1}x_{42}}\rightarrow 0}}{{x_{33}\rightarrow 1},{x_{44}\rightarrow 1}}}}{{\overset{*}{=}\begin{pmatrix}1 & {- 1} & 0 & {- 1} \\{- 1} & 1 & 1 & 0 \\0 & 1 & 1 & 1 \\{- 1} & 0 & 1 & 1\end{pmatrix}};{{\underset{4}{a}}_{k\overset{\_}{l}}{_{\underset{{x_{43}\rightarrow 0},{y_{43}\rightarrow 0}}{\underset{{x_{41}\rightarrow{- 1}},{x_{42}\rightarrow 0}}{{x_{33}\rightarrow 1},{x_{44}\rightarrow 1}}}}{\overset{*}{=}{\begin{pmatrix}{- 1} & {- 1} & {- 1} & {- 1} \\{- 1} & 1 & 0 & 0 \\{- 1} & 0 & 1 & 0 \\{- 1} & 0 & 0 & 1\end{pmatrix}.}}}}}}}}}}}}} & (71)\end{matrix}$

Once again we pick a solution and then illustrate a few of the basicquantities that are derived for the DFT in Schouten's notation. Forgenerality, consider the complex-valued solution of (71) for the a_(kl):

$\begin{matrix}{{{a_{\overset{\_}{k}l} \equiv {\underset{2}{a}}_{\overset{\_}{k}l}}\overset{*}{=}\begin{pmatrix}1 & {- i} & {2i} & {- i} \\i & 1 & i & 0 \\{{- 2}i} & {- i} & 1 & {- i} \\i & 0 & i & 1\end{pmatrix}},} & (72)\end{matrix}$

and then the inverse is given by

$\begin{matrix}{{\overset{- 1}{a}}_{\overset{\_}{k}l} = {a^{k\overset{\_}{l}}\overset{*}{=}{\frac{1}{7}{\begin{pmatrix}1 & {{- 2} - i} & {{- 2} + {2i}} & {{- 2} - i} \\{{- 2} + {2i}} & 5 & {2 + i} & {- 2} \\{{- 2} - {2i}} & {2 - i} & 1 & {2 - i} \\{{- 2} + i} & {- 2} & {2 + i} & 5\end{pmatrix}.}}}} & (73)\end{matrix}$

To illustrate some specific quantities we define a general complex“data” vector:

f ^(n)=(z ¹ ,z ² ,z ³ ,z ⁴),  (74)

and then the DFT transformation (9) is given by

$\begin{matrix}{F^{m} = {{W_{\cdot \; n}^{m}f^{n}} = {\frac{1}{2}{\begin{pmatrix}{z^{1} + z^{2} + z^{3} + z^{4}} \\{z^{1} - {iz}^{2} - z^{3} + {iz}^{4}} \\{z^{1} - z^{2} + z^{3} - z^{4}} \\{z^{1} + {iz}^{2} - z^{3} - {iz}^{4}}\end{pmatrix}.}}}} & (75)\end{matrix}$

Applying the inverse DFT transformation shows that as expected:

$\begin{matrix}\begin{matrix}{{{\overset{\sim}{W}}_{\cdot m}^{n}F^{m}} = {\frac{1}{4}\begin{pmatrix}1 & 1 & 1 & 1 \\1 & i & {- 1} & {- i} \\1 & {- 1} & 1 & {- 1} \\1 & {- i} & {- 1} & i\end{pmatrix}\begin{pmatrix}{z^{1} + z^{2} + z^{3} + z^{4}} \\{z^{1} - {iz}^{2} - z^{3} + {iz}^{4}} \\{z^{1} - z^{2} + z^{3} - z^{4}} \\{z^{1} + {iz}^{2} - z^{3} - {iz}^{4}}\end{pmatrix}}} \\{= \left( {z^{1},z^{2},z^{3},z^{4}} \right)} \\{= {f^{n}.}}\end{matrix} & (76)\end{matrix}$

Other quantities are derived from the f^(n) using the a _(kl) of (72),e.g., consider

$\begin{matrix}{{f_{\overset{\_}{k}} = {{a_{\overset{\_}{k}n}f^{n}} = {\frac{1}{2}\begin{pmatrix}{z_{1} - {iz}_{2} + {2\; {iz}_{3}} - {iz}_{4}} \\{{iz}_{1} + z_{2} + {iz}_{3}} \\{{{- 2}\; {iz}_{1}} - {iz}_{2} + z_{3} - {iz}_{4}} \\{{iz}_{1} + {iz}_{3} + z_{4}}\end{pmatrix}}}},} & (77)\end{matrix}$

and we further check that f _(n) =

, by calculating:

$\begin{matrix}\begin{matrix}{{a_{\overset{\_}{k}n}\overset{\_}{\left( f^{k} \right)}} = {a_{\overset{\_}{k}n}{\overset{\_}{f}}^{\overset{\_}{k}}}} \\{= {\overset{\_}{f}}_{n}} \\{= {\frac{1}{2}{\begin{pmatrix}{{\overset{\_}{z}}_{1} + {i{\overset{\_}{z}}_{2}} - {2\; i{\overset{\_}{z}}_{3}} + {i{\overset{\_}{z}}_{4}}} \\{{{- i}{\overset{\_}{z}}_{1}} + {\overset{\_}{z}}_{2} - {i{\overset{\_}{z}}_{3}}} \\{{{- 2}\; i{\overset{\_}{z}}_{1}} + {i{\overset{\_}{z}}_{2}} + {\overset{\_}{z}}_{3} + {i{\overset{\_}{z}}_{4}}} \\{{{- i}{\overset{\_}{z}}_{1}} - {i{\overset{\_}{z}}_{3}} + {\overset{\_}{z}}_{4}}\end{pmatrix}.}}}\end{matrix} & (78)\end{matrix}$

Comparing f _(n) in (77) to the complex conjugation of f _(n) in (78) wesee that as expected:

f _(n) =

,   (79)

demonstrating consistency using the hybrid-index notation.Correspondingly, we can lower indices on the F^(m) to get

$\begin{matrix}\begin{matrix}{F_{\overset{\_}{k}} = {a_{\overset{\_}{k}m}F^{m}}} \\{= {\frac{1}{2}{\begin{pmatrix}{z^{1} + {\left( {1 - {2i}} \right)z^{2}} + {\left( {1 + {4i}} \right)z^{3}} + {\left( {1 - {2i}} \right)z^{4}}} \\{{\left( {1 + {2i}} \right)z^{1}} - {iz}^{2} - {\left( {1 - {2i}} \right){iz}^{3}} + {iz}^{4}} \\{{\left( {1 - {4i}} \right)z^{1}} - {\left( {1 + {2i}} \right)z^{2}} + z^{3} - {\left( {1 + {2i}} \right)z^{4}}} \\{{\left( {1 + {2i}} \right)z^{1}} + {iz}^{2} + {\left( {{2i} - 1} \right)z^{3}} - {iz}^{4}}\end{pmatrix}.}}}\end{matrix} & (80)\end{matrix}$

Parseval's theorem is invariant as expected:

$\begin{matrix}{{{{\underset{2}{a}}_{\overset{\_}{m}n}{\overset{\_}{F}}^{\overset{\_}{m}}F^{n}} = {{{\underset{2}{a}}_{\overset{\_}{k}l}{\overset{\_}{f}}^{\overset{\_}{k}}f^{l}} = {{z^{1}}^{2} + {z^{2}}^{2} + {z^{3}}^{2} + {z^{4}}^{2}}}},} & (81)\end{matrix}$

but the Euclidean form in (81) deriving from this particular

${\underset{2}{a}}_{\overset{\_}{m}n}$

is not necessarily expected. As pointed out earlier, Parseval's theoremcan take different forms for the different a _(kl) solutions. E.g.,considering solution 4 of the a _(kl) in (71) we have

$\begin{matrix}\begin{matrix}{{{\underset{4}{a}}_{\overset{\_}{m}n}{\overset{\_}{F}}^{\overset{\_}{m}}F^{n}} = {{\underset{4}{a}}_{\overset{\_}{k}l}{\overset{\_}{f}}^{\overset{\_}{k}}f^{l}}} \\{{= {{- {z^{1}}^{2}} + {z^{2}}^{2} + {z^{3}}^{2} + {z^{4}}^{2} - {2{z^{1}\left( {z^{2} + z^{3} + z^{4}} \right)}}}},}\end{matrix} & (82)\end{matrix}$

which is invariant but not real valued.

Eigenvalue-Eigenvector Problem

Given the formalism above, we can examine the eigenvalue-eigenvector(eigen-problem) from a more general viewpoint than what is normallydiscussed. The additional generality stems from the possibility that themetric is not necessarily the Kronecker delta, constrained only by itsinvariance under unitary transformations and hermitian symmetry. We haveshown from a general viewpoint that the a _(kl)≐δ _(kl) solution is animportant “base” solution without additional worries that perhapssomething was overlooked. We have also shown that other solutions forthe a _(kl) exist based simply on unitary invariance, and thus, moregeneral calculations of the DFT and its various quantities are easilyextended to arbitrary basis sets. Some advantages in calculation may bepossible in these other general basis sets and we will reserve thistopic for a future set of notes. For now, we examine a specific basisset that follows by solution to the eigen-problem as discussed below.

Before examining the specifics of the eigen-problem for the DFT, somegeneralities are discussed. The eigenvalue-eigenvector problem seekssolutions to the following set equations:

$\begin{matrix}{{{T_{\overset{\_}{m}n}{\underset{a}{v}}^{n}} = {\underset{a}{\sigma}{\underset{a}{v}}_{\overset{\_}{m}}}},} & (83)\end{matrix}$

for the vectors

${\underset{a}{v}}^{n}$

of a general rank-2 quantity, T _(mn). The scalar values,

$\underset{a}{\sigma},$

are the eigenvalues and the

${\underset{a}{v}}^{n}$

are the eigenvectors. The “a” index below the kernel object, in thiscase

$,\underset{a}{v}$

and

$\underset{a}{\sigma}$

does not refer to a component, but rather, the a-th enigenvector,

$\underset{a}{v}$

corresponding to the a-th eigenvalue,

$\underset{a}{\sigma}$

as separate and distinct solutions. We therefore choose indices from theset, {a, b, c, d}, to label the kernel objects corresponding to distincteigen-solutions, and then reserve the index range, {k, l, m, n}, tolabel components of a given basis.

Continuing to the DFT application we substitute W _(mn) of (26) into(83):

$\begin{matrix}{{{W_{\overset{\_}{m}n}{\underset{a}{v}}^{n}} = {\underset{a}{\sigma}{\underset{a}{v}}_{\overset{\_}{m}}}},} & (84)\end{matrix}$

and make special note that we have switched to a common kernel symbolfor the Fourier transform pair, (F^(m), f^(n))→(v _(m) , v^(n)), withthe substitutions:

v ^(n) ≡f ^(n) and F _(m) ≡v _(m) .   (85)

Looking ahead, we mention now that the solution to (84) results in atransformation of the W _(mn) into a diagonal basis by use of thesimilarity transform. In the context of the present application, andnoting expression (9) for the Fourier transform, (84) then has aninteresting geometrical interpretation as solving for the Fouriertransform pair, (v _(m) , v^(n)), that are “parallel” to one another, orequivalently, left invariant under the action of the DFT (aside from thescaling factor σ). To give greater emphasis to this last comment,consider the overall structure of the DFT transformation in a diagonalbasis for the 4×4 case (for now using the notation of (9) to emphasizethis last point):

$\begin{matrix}\begin{matrix}{{W_{\overset{\_}{m}n}{\underset{a}{v}}^{n}} = \left. {\underset{a}{\sigma}{\underset{a}{v}}_{\overset{\_}{m}}}\Rightarrow{\begin{pmatrix}\underset{1}{\sigma} & 0 & 0 & 0 \\0 & \underset{2}{\sigma} & 0 & 0 \\0 & 0 & \underset{3}{\sigma} & 0 \\0 & 0 & 0 & \underset{4}{\sigma}\end{pmatrix}\begin{pmatrix}f^{1^{\prime}} \\f^{2^{\prime}} \\f^{3^{\prime}} \\f^{4^{\prime}}\end{pmatrix}} \right.} \\{= \begin{pmatrix}{\underset{1}{\sigma}f^{1^{\prime}}} \\{\underset{2}{\sigma}f^{2^{\prime}}} \\{\underset{3}{\sigma}f^{2^{\prime}}} \\{\underset{4}{\sigma}f^{2^{\prime}}}\end{pmatrix}} \\{{= {\underset{a}{\sigma}\begin{pmatrix}F_{1^{\prime}} \\F_{2^{\prime}} \\F_{3^{\prime}} \\F_{4^{\prime}}\end{pmatrix}}},}\end{matrix} & (86)\end{matrix}$

where the primes on indices denote the data and its transform in thediagonal basis as discussed further below. Thus as mentioned above, theFourier transform pair are co-linear when calculated in the diagonalbasis:

$\begin{matrix}{{\begin{pmatrix}F_{1^{\prime}} \\F_{2^{\prime}} \\F_{3^{\prime}} \\F_{4^{\prime}}\end{pmatrix} = \begin{pmatrix}{\left( {\underset{1}{\sigma}/\underset{a}{\sigma}} \right)f^{1^{\prime}}} \\{\left( {\underset{2}{\sigma}/\underset{a}{\sigma}} \right)f^{2^{\prime}}} \\{\left( {\underset{3}{\sigma}/\underset{a}{\sigma}} \right)f^{3^{\prime}}} \\{\left( {\underset{4}{\sigma}/\underset{a}{\sigma}} \right)f^{4^{\prime}}}\end{pmatrix}},} & (87)\end{matrix}$

differing only by a scalar factor. This example illustrates thegeometrical significance of the DFT transform in a diagonal basis andmotivates using the same kernel symbol for the Fourier transform pair in(85), since the quantities are co-linear. This difference in notation isto be contrasted with the approach taken so far where the kernel symbolshave been kept distinct, and the W_(·n) ^(m) transformation was regardedas a point transformation.

Continuing with the eigen-problem calculation we use the metric tensorto show that (84) is equivalently expressed in the form:

$\begin{matrix}{{{\left( {W_{\overset{\_}{m}n} - {\underset{a}{\sigma}a_{\overset{\_}{m}n}}} \right){\underset{a}{v}}^{n}} = 0},} & (88)\end{matrix}$

or also as

$\begin{matrix}{{\left( {W_{\cdot n}^{m} - {\underset{a}{\sigma}\delta_{\cdot n}^{m}}} \right){\underset{a}{v}}^{n}} = 0.} & (89)\end{matrix}$

And thus the eigen-problem, when cast in the form of (89), does notdepend on a specific choice for the metric. Thus we solve the DFTeigen-problem in the form of (89). The eigenvalues for the W_(·n) ^(m)are then solutions to

$\begin{matrix}{{\det \left( {W_{\cdot n}^{m} - {\underset{a}{\sigma}\delta_{\cdot n}^{m}}} \right)} = 0.} & (90)\end{matrix}$

Eigen-Problem for the DFT Matrix

The eigenvector-eigenvalue problem for the DFT matrix in the form of(89) has been examined extensively and so we review some of these basicresults before examining the diagonal DFT calculation in detail. Since

W⁴=I,  (91)

the eigenvalues for the DFT matrix for N>3 are non-unique and aremultiples of the so-called “roots of unity:”

$\begin{matrix}{{{\underset{a}{\sigma} = ^{{\pi}\; {a/2}}};{a = 1}},\ldots \mspace{14mu},{\left. {4;}\Rightarrow\underset{a}{\sigma} \right. \in {\left\{ {{\pm 1},{\pm i}} \right\}.}}} & (92)\end{matrix}$

Specifically, (92) then becomes

$\begin{matrix}{{\underset{1}{\sigma} = 1},{\underset{2}{\sigma} = {- 1}},{\underset{3}{\sigma} = i},{\underset{4}{\sigma} = {- {i.}}}} & (93)\end{matrix}$

The eigenvalues of the DFT transform are not all distinct for N≧4, sothere are many possible sets of eigenvectors satisfying (89). As aconsequence of having these non-unique eigenvalues, the correspondingeigenvectors for these degenerate eigenvalues are not necessarilyorthogonal. The importance of having an orthogonal set of eigenvectorsfor the present application is not known, but should this turn out to bean issue, the Gramm-Schmidt procedure can always be applied or anorthogonal set can be obtained. Further consequences and examples ofthis analysis for higher dimensional DFTs will be discussed in a futureset of notes. For now we examine a low-dimensional example to map outsome of the basic features.

Eigenvalues

Considering M=N=4, the DFT transform and its inverse is given in (68)which is reproduced below:

${W_{\cdot k}^{m} = {\frac{1}{2}\begin{pmatrix}1 & 1 & 1 & 1 \\1 & {- i} & {- 1} & i \\1 & {- 1} & 1 & {- 1} \\1 & i & {- 1} & {- i}\end{pmatrix}}};$${\overset{- 1}{W}}_{\cdot k}^{m} = {{W\mspace{14mu} \%_{\cdot k}^{m}} = {{\overset{\_}{W}}_{k}^{\cdot m} = {\frac{1}{2}{\begin{pmatrix}1 & 1 & 1 & 1 \\1 & i & {- 1} & {- i} \\1 & {- 1} & 1 & {- 1} \\1 & {- i} & {- 1} & i\end{pmatrix}.}}}}$

Considering solutions of (90) we have

${{\det \left( {W_{\cdot n}^{m} - {\underset{a}{\sigma}\delta_{\cdot n}^{m}}} \right)} = 0},$

and for the 4×4 example has the solution:

$\begin{matrix}{{{\left( {\underset{a}{\sigma} - 1} \right)^{2}\left( {\underset{a}{\sigma} + 1} \right)\left( {\underset{a}{\sigma} + i} \right)} = \left. 0\Rightarrow\left\{ {{\underset{1}{\sigma} = 1},{\underset{2}{\sigma} = {- 1}},{\underset{3}{\sigma} = 1},{\underset{4}{\sigma} = {- i}}} \right\} \right.},} & (94)\end{matrix}$

with one degenerate eigenvalue,

$\underset{1}{\sigma} = {\underset{3}{\sigma} = 1.}$

Multiplicity of Eigenvalues for the DFT Matrix

As discussed above, the significance of (91) for the DFT eigenvalueproblem is that degenerate eigenvalues appear for N≧4. The“multiplicities” of these eigenvalues are defined as the number of timeseach eigenvalue appears. For the 4×4 example in (94) we therefore have

m ₁=2,m _(2′)=1,m ₃=1,m ₄=0,   (95)

where the m_(k) ordering corresponds to the ordering of (93). Theseresults can be generalized to higher dimensions by the following Tableof multiplicities for each eigenvalue:

TABLE 1 M = N m₁${{for}\mspace{14mu} \underset{1}{\sigma}} = {1\text{:}}$ m₂${{for}\mspace{14mu} \underset{2}{\sigma}} = {{- 1}\text{:}}$ m₃${{for}\mspace{14mu} \underset{3}{\sigma}} = {i\text{:}}$ m₄${{for}\mspace{14mu} \underset{4}{\sigma}} = {{- i}\text{:}}$ 4m m +1 m m m − 1 4m + 1 m + 1 m m m 4m + 2 m + 1 m + 1 m m 4m + 3 m + 1 m + 1m + 1 m

By inspection of Table 1, the 4×4 case corresponds to the m=1 value inthe top row giving the multiplicities of eigenvalues listed in (95).Similarly, it can be shown that for the 5×5 case, the eigenvalues aresolutions of

$\begin{matrix}{{{\left( {\underset{a}{\sigma} - 1} \right)^{2}\left( {\underset{a}{\sigma} + 1} \right)\left( {\underset{a}{\sigma^{2}} + 1} \right)} = \left. 0\Rightarrow\left\{ {{\underset{1}{\sigma} = 1},{\underset{2}{\sigma} = {- 1}},{\underset{3}{\sigma} = i},{\underset{4}{\sigma} = {- i}},{\underset{5}{\sigma} = 1}} \right\} \right.},} & (96)\end{matrix}$

corresponding to the 2^(nd) row of the Table with m=1. From (96) (orequivalently Table 1) the multiplicities for the 5×5 case are thus

m ₁=2,m ₂=1,m ₃=1,m ₄=1.  (97)

To illustrate a higher dimensional example, consider say M=N=131=4.32+3,corresponding to the third row of the Table giving the followingmultiplicities:

m ₁=33,m ₂=22,m ₃=22,m ₄=32.   (98)

Eigenvectors

Given the eigenvalues in (94) the eigenvectors corresponding to eacheigenvalue are constructed using (89). It is instructive to considerthis process for the eigenvectors corresponding to the degenerate firstand third eigenvalues,

$\underset{1}{\sigma} = {\underset{3}{\sigma} = 1.}$

The first eigenvector is obtained as a solution to the set of linearequations:

$\begin{matrix}{{{{\left( {W_{\cdot n}^{m} - {\underset{a}{\sigma}\delta_{\cdot n}^{m}}} \right){\underset{a}{v}}^{n}}_{a\rightarrow 1}} = {\left. 0\Rightarrow{\begin{pmatrix}{1 - \underset{1}{\sigma}} & 1 & 1 & 1 \\1 & {{- i} - \underset{1}{\sigma}} & {- 1} & i \\1 & {- 1} & {1 - \underset{1}{\sigma}} & {- 1} \\1 & i & {- 1} & {{- i} - \underset{1}{\sigma}}\end{pmatrix}\begin{pmatrix}{\underset{1}{v}}^{1} \\{\underset{1}{v}}^{2} \\{\underset{1}{v}}^{3} \\{\underset{1}{v}}^{4}\end{pmatrix}} \right. = \begin{pmatrix}0 \\0 \\0 \\0\end{pmatrix}}},} & (99)\end{matrix}$

which simplifies to

$\begin{matrix}{\begin{pmatrix}{{- {\underset{1}{v}}^{1}} + {\underset{1}{v}}^{2} + {\underset{1}{v}}^{3} + {\underset{1}{v}}^{4}} \\{{\underset{1}{v}}^{1} - {\left( {2 + i} \right){\underset{1}{v}}^{2}} - {{\underset{1}{v}}^{3}i{\underset{1}{v}}^{4}}} \\{{\underset{1}{v}}^{1} - {\underset{1}{v}}^{2} - {\underset{1}{v}}^{3} - {\underset{1}{v}}^{4}} \\{{\underset{1}{v}}^{1} + {i{\underset{1}{v}}^{2}} - {\underset{1}{v}}^{3} - {\left( {2 + i} \right){\underset{1}{v}}^{4}}}\end{pmatrix} = {\begin{pmatrix}0 \\0 \\0 \\0\end{pmatrix}.}} & (100)\end{matrix}$

The solution to (100) for

${\underset{1}{v}}^{n}$

is therefore given in terms of two arbitrary components, say

${\underset{1}{v}}^{3}\mspace{14mu} {and}\mspace{14mu} {\underset{1}{v}}^{4}\text{:}$

$\begin{matrix}{\left( {{\underset{1}{v}}^{1},{\underset{1}{v}}^{2},{\underset{1}{v}}^{3},{\underset{1}{v}}^{4}} \right) = {\left( {{{\underset{1}{v}}^{3} + {2{\underset{1}{v}}^{4}}},{\underset{1}{v}}^{4},{\underset{1}{v}}^{3},{\underset{1}{v}}^{4}} \right).}} & (101)\end{matrix}$

We point out that two components are required because of the two-folddegeneracy of the

$\underset{1}{\sigma} = {\underset{3}{\sigma} = 1}$

eigenvalue. In general, the eigenvectors corresponding to uniqueeigenvalues are specified in terms of a single component. Continuing wearbitrarily choose

${\underset{1}{v}}^{4} = 1$

and

${{\underset{1}{v}}^{3} = 0},$

to obtain a realization for the first eigenvector:

$\begin{matrix}\begin{matrix}{\left( {{\underset{1}{v}}^{1},{\underset{1}{v}}^{2},{\underset{1}{v}}^{3},{\underset{1}{v}}^{4}} \right) = {\left( {{{\underset{1}{v}}^{3} + {2{\underset{1}{v}}^{4}}},{\underset{1}{v}}^{4},{\underset{1}{v}}^{3},{\underset{1}{v}}^{4}} \right)_{\underset{{\underset{1}{v}}^{3}\rightarrow 0}{{\underset{1}{v}}^{4}\rightarrow 1}}\left. \Rightarrow\underset{1}{v^{n}} \right.}} \\{= {\left( {2,1,0,1} \right).}}\end{matrix} & (102)\end{matrix}$

Note also that since,

${\underset{1}{\sigma} = {\underset{3}{\sigma} = 1}},$

are degenerate, the solution to (99) is trivially the same as (100).Therefore, a solution for

${\underset{3}{v}}^{n}$

follows immediately by choosing complementary but arbitrary numericalvalues for the

${\underset{1}{v}}^{3}$

and

${\underset{1}{v}}^{4}$

in (102), e.g.,

${\underset{3}{v}}^{4} = 0$

and

$\begin{matrix}{{\underset{3}{v}}^{3} = {{1:\left( {{\underset{3}{v}}^{1},{\underset{3}{v}}^{2},{\underset{3}{v}}^{3},{\underset{3}{v}}^{4}} \right)} = {{\left( {{{\underset{3}{v}}^{3} + {2{\underset{3}{v}}^{4}}},{\underset{3}{v}}^{4},{\underset{3}{v}}^{3},{\underset{3}{v}}^{4}} \right)_{{\underset{1}{v}}^{3}\rightarrow 1}^{\text{?}}\left. \Rightarrow{\underset{3}{v}}^{n} \right.} = {{\left( {1,0,1,0} \right).\text{?}}\text{indicates text missing or illegible when filed}}}}} & (103)\end{matrix}$

Similarly, we can solve for the remaining eigenvectors corresponding totheir eigenvalues,

$\underset{2}{\sigma} = {- 1}$

and

$\underset{4}{\sigma} = {- {i.}}$

In summary, a set of engenvectors (but not the only set) for the 4×4 DFTtransformation is summarized below:

$\begin{matrix}{{{\underset{1}{v}}^{n} = \begin{pmatrix}2 \\1 \\0 \\1\end{pmatrix}},{{\underset{2}{v}}^{n} = \begin{pmatrix}{- 1} \\1 \\1 \\1\end{pmatrix}},{{\underset{3}{v}}^{n} = \begin{pmatrix}1 \\0 \\1 \\0\end{pmatrix}},{{\underset{4}{v}}^{n} = {\begin{pmatrix}0 \\{- 1} \\0 \\1\end{pmatrix}.}}} & (104)\end{matrix}$

Since two of the eigenvalues are degenerate, (102) and (104)collectively do not form an orthogonal system as can be verified byforming the inner product between any two vectors using the “base”solution,

${\underset{1}{a}}_{\overset{\_}{k}l}\overset{*}{=}\delta_{\overset{\_}{k}l}$

of (71):

$\begin{matrix}{{{\langle{\underset{m}{v}\underset{n}{v}}\rangle} = {{{\underset{1}{a}}_{\overset{\_}{k}l}{\overset{\_}{\underset{m}{v}}}^{\overset{\_}{k}}{\underset{n}{v}}^{l}} = {{\delta_{\overset{\_}{k}l}{\overset{\_}{\underset{m}{v}}}^{\overset{\_}{k}}{\underset{n}{v}}^{l}} = \begin{pmatrix}3 & 0 & 1 & 0 \\0 & 2 & 0 & 0 \\1 & 0 & 1 & 0 \\0 & 0 & 0 & 1\end{pmatrix}}}},} & (105)\end{matrix}$

where the off-diagonal entries correspond to the

$\underset{1}{\sigma} = {\underset{3}{\sigma} = 1}$

value. For reference, we can still construct an orthonormal set from(104):

$\begin{matrix}{{{\underset{1}{u}}^{n} = {\frac{1}{\sqrt{6}}\begin{pmatrix}2 \\1 \\0 \\1\end{pmatrix}}},{{\underset{2}{u}}^{n} = {{- \frac{1}{2}}\begin{pmatrix}{- 1} \\1 \\1 \\1\end{pmatrix}}},{{\underset{3}{u}}^{n} = {\frac{1}{2\sqrt{3}}\begin{pmatrix}1 \\{- 1} \\3 \\{- 1}\end{pmatrix}}},{{\underset{4}{u}}^{n} = {\frac{1}{\sqrt{2}}{\begin{pmatrix}0 \\{- 1} \\0 \\1\end{pmatrix}.}}}} & (106)\end{matrix}$

and verify as expected that

$\begin{matrix}{{\delta_{\overset{\_}{k}l}{\overset{\_}{\underset{m}{u}}}^{\overset{\_}{k}}{\underset{n}{u}}^{l}} = {\delta_{mn}.}} & (107)\end{matrix}$

Similarity Transform

Given the eigenvectors of (104) we form the “similarity” transform usingthe eigenvectors as columns. We denote this as a transformation of basisand we find by construction that the similarity transform, S_(m) ^(m)′,and its inverse,

${\overset{- 1}{S}}_{m}^{m^{\prime}},$

are given by

$\begin{matrix}{{S_{m}^{m^{\prime}} = \begin{pmatrix}2 & {- 1} & 1 & 0 \\1 & 1 & 0 & {- 1} \\0 & 1 & 1 & 0 \\1 & 1 & 0 & 1\end{pmatrix}};{{\overset{- 1}{S}}_{m}^{m^{\prime}} = {\frac{1}{4}{\begin{pmatrix}1 & 1 & {- 1} & 1 \\{- 1} & 1 & 1 & 1 \\1 & {- 1} & 3 & {- 1} \\0 & {- 2} & 0 & 2\end{pmatrix}.}}}} & (108)\end{matrix}$

When using the kernel-index method for transformations of basis, sayx^(m)′=A_(m) ^(m)′x^(m). that the A_(m) ^(m)′ quantities do not use dotsabove or below an index to mark index placement. This differs from thenotation for point transformations, which does use the dots e.g.,F^(m)=W_(·n) ^(m)f^(n). The reason for this is to discourage the“raising” or “lowering” of indices on the A_(m) ^(m)′ to avoid aninconsistency when using the metric tensor. This inconsistency can occurbecause the metric a _(kl) also transforms under the A_(m) ^(m)′, andthus has different components in the new coordinate system. For pointtransformations however, no coordinate transformation occurs and thusindices can be raised or lowered using a single (unambiguous) metric.

We then verify that applying the similarity transformation to the W_(·n)^(m) gives as expected:

$\begin{matrix}{{{SWS}^{- 1} \equiv {S_{m}^{m^{\prime}}W_{\cdot n}^{m}{\overset{- 1}{S}}_{n^{\prime}}^{n}}} = {W_{\cdot n^{\prime}}^{m^{\prime}} = {\begin{pmatrix}\underset{1}{\sigma} & 0 & 0 & 0 \\0 & \underset{2}{\sigma} & 0 & 0 \\0 & 0 & \underset{3}{\sigma} & 0 \\0 & 0 & 0 & \underset{4}{\sigma}\end{pmatrix} = {\begin{pmatrix}1 & 0 & 0 & 0 \\0 & {- 1} & 0 & 0 \\0 & 0 & 1 & 0 \\0 & 0 & 0 & {- i}\end{pmatrix}.}}}} & (109)\end{matrix}$

As a check on index ordering we run the following line of Mathematicacode with Si defined as the inverse of the Sin (108):

Table [Sum[S[[m,ump]]W[[m,n]]Si[[lnp,n]], {n,l,M}, {m,l,M}], {ump,l,M},{lnp,l,M}]  (110)

where ump denotes “upper-m-prime” and lnp is “lower-m-prime.”

It is of further interest to point out that if the similarity transformis constructed using the orthonormal eigenvectors of (106) as columns(and correspondingly it inverse), then the application of the“orthonormal” similarity transform to the W_(·n) ^(m) also gives (109).Therefore, the similarity transform that leaves the W_(·n) ^(m) indiagonal form is certainly not unique, so many possible choices can bemade for the S_(n) ^(m)′. So as indicated by calculations thus far, itdoes not seem to matter whether or not the eigenvectors forming theS_(n) ^(m)′ are orthogonal, and this can greatly simplify the initialcalculations needed for calculating in the diagonal basis.

DFT Calculation in Diagonal Form

The diagonal form of the DFT in (109) shows that the DFT calculation canbe performed both efficiently and elegantly as simply a multiplicationof the entries on the main diagonal (the eigenvalues), to thecorresponding data components that have been pre-multiplied using thesimilarity matrix. This means that after pre-multiplying the data, theDFT transform is reduced to a vector product consisting of N operationsrather than a matrix multiply with N² operations (or N log N when usingthe FFT). In reality though, this reduction in computation cost is onlypartly true in terms of the performance gains, i.e., the diagonal DFTrepresentation does obviously lead to a speedy calculation for one partof the DFT operation, but in general, this gain is offset by the addedcomputational cost of having to pre-transform the data and thenback-transform to bring the Fourier transform into the original basis asillustrated below. But the key is that these latter operations may notneed to be performed at every step in an iterative transformapplication, such as with phase retrieval. Additional details on thiscalculation approach are outlined below.

To begin, we write down the DFT calculation in the diagonal basis afterswitching the notation back to distinct kernel symbols for the data andits transform:

f ^(m) ′=W _(·n) ^(m) ′f ^(n)′,   (111)

where the m′ represents the transformation of basis under the action ofthe similarity transform as illustrated above in (109). Since the onlynon-zero values for the W_(·n′) ^(m)′ are along the diagonal, we form avector,

$\underset{m^{\prime}}{\sigma},$

form the eigenvalues (i.e., the diagonal entries of the W_(·n′) ^(m)′)as in (94). Continuing with the calculation of (111, the starting data,f_(k), are pre-multiplied using the inverse similarity transform,

${\overset{- 1}{S}}_{n}^{n^{\prime}},$

to give

$\begin{matrix}{f^{n^{\prime}} = {{\overset{- 1}{S}}_{n}^{n^{\prime}}{f^{n}.}}} & (112)\end{matrix}$

The diagonal DFT is then the product of the N components of the

$\underset{m^{\prime}}{\sigma}$

and the f^(n)′ defined by (112):

$\begin{matrix}{{F^{m^{\prime}} = {\underset{m^{\prime}}{\sigma} \cdot f_{m^{\prime}}}},} & (113)\end{matrix}$

where no sum is implied over the (m′,m′) pair. To complete thecalculation, the “diagonal” DFT components in (113) need to betransformed back to the original DFT basis, F^(m). To obtain this resultwe express the LHS of (111) in terms of components in the originalbasis:

$\begin{matrix}{{F^{m^{\prime}} = {{\overset{- 1}{S}}_{m}^{m^{\prime}}F^{m}}},} & (114)\end{matrix}$

and therefore the original Fourier transform components are readilyobtained by back-transforming the F^(m)′ of (114) using the similaritytransform:

$\begin{matrix}{F^{n} = {{S_{n^{\prime}}^{n}F^{n^{\prime}}} = {{S_{n^{\prime}}^{n}{\overset{- 1}{S}}_{m}^{n^{\prime}}F^{m}} = {\delta_{m}^{n}{F^{m}.}}}}} & (115)\end{matrix}$

The inverse DFT is also calculated in a diagonal basis since the inverseof a diagonal matrix is found by populating its diagonal entries withthe inverse of the individual components. Therefore, given the W_(·n)^(m)′,the inverse is simply (no sum is implied on the RHS over the (m′,m′) pair):

$\begin{matrix}{{{\overset{- 1}{W}}_{\bullet \; n^{\prime}}^{m^{\prime}} = {\left( {1/\underset{m^{\prime}}{\sigma}} \right)\delta_{n^{\prime}}^{m^{\prime}}}},} & (116)\end{matrix}$

and then inverse DFT transformation in the diagonal basis gives asexpected:

$\begin{matrix}{{{\left( {1/\underset{m^{\prime}}{\sigma}} \right) \cdot F^{m^{\prime}}} = {{\left( {1/\underset{m^{\prime}}{\sigma}} \right) \cdot \underset{m^{\prime}}{\sigma} \cdot f^{m^{\prime}}} = f^{m^{\prime}}}},} & (117)\end{matrix}$

The overall point is that for applications where it is feasible toperform the DFT in a diagonal basis, the similarity transform and itsinverse can be pre-calculated. The multiplication of the similaritytransform with the data can also be pre-calculated. Then at anappropriate point, the DFT results can be transformed back to theoriginal basis using the inverse similarity transform, preferably at theend of an iterative loop. The approach is summarized in FIG. 2.

Example: 4×4 Diagonal DFT Calculation

The steps above are illustrated for the 4×4 example. Consider thegeneral complex data vector of (74):

f ^(n)=(z ¹ ,z ² ,z ³ ,z ⁴).  (118)

The inverse similarity transform of (108) times the data vector (118) isgiven by

$\begin{matrix}{{f^{m^{\prime}} = {{{\overset{- 1}{S}}_{m}^{m^{\prime}}f^{m}} = {\frac{1}{4}\begin{pmatrix}{z_{1} + z_{2} - z_{3} + z_{4}} \\{{- z_{1}} + z_{2} + z_{3} + z_{4}} \\{z_{1} - z_{2} + {3z_{3}} - z_{4}} \\{2\left( {{- z_{2}} + z_{4}} \right)}\end{pmatrix}}}},} & (119)\end{matrix}$

and then from (94):

$\begin{matrix}{\underset{m^{\prime}}{\sigma} = {\left( {1,{- 1},1,{- i}} \right).}} & (120)\end{matrix}$

The DFT in the transformed basis follows by multiplication of theeigenvalues in (120) with the f^(m)′ components in (119) to give theFourier transform in the diagonal basis:

$\begin{matrix}{F^{m^{\prime}} = {{\underset{m^{\prime}}{\sigma} \cdot f^{m}} = {\frac{1}{4}{\begin{pmatrix}{z_{1} + z_{2} - z_{3} + z_{4}} \\{z_{1} - z_{2} - z_{3} - z_{4}} \\{z_{1} - z_{2} + {3z_{3}} - z_{4}} \\{{- 2}{i\left( {{- z_{2}} + z_{4}} \right)}}\end{pmatrix}.}}}} & (121)\end{matrix}$

Finally, we back-transform the F^(m)′ components in (121) by multiplyingwith the similarity transformation of (108) to give the DFT result inthe original non-diagonal basis:

$\begin{matrix}\begin{matrix}{F^{m} = {S_{m^{\prime}}^{m}F^{m^{\prime}}}} \\{= {\begin{pmatrix}2 & {- 1} & 1 & 0 \\1 & 1 & 0 & {- 1} \\0 & 1 & 1 & 0 \\1 & 1 & 0 & 1\end{pmatrix}\frac{1}{4}\begin{pmatrix}{z_{1} + z_{2} - z_{3} + z_{4}} \\{z_{1} - z_{2} - z_{3} - z_{4}} \\{z_{1} - z_{2} + {3z_{3}} - z_{4}} \\{{- 2}{i\left( {{- z_{2}} + z_{4}} \right)}}\end{pmatrix}}} \\{{= {\frac{1}{2}\begin{pmatrix}{z^{1} + z^{2} + z^{3} + z^{4}} \\{z^{1} - {iz}^{2\;} - z^{3} + {iz}^{4}} \\{z^{1} - z^{2} + z^{3} - z^{4}} \\{z^{1} + {iz}^{2} - z^{3} - {iz}^{4}}\end{pmatrix}}},}\end{matrix} & (122)\end{matrix}$

in agreement with the original DFT transform, F^(m)=W_(·n) ^(m) f^(n) of(75). Similarly, the inverse transformation can be calculated in thediagonal basis as discussed above using (117).

The diagonal form of the DFT in (113) shows that the calculation can beperformed, efficiently and elegantly as simply a vector product. Thismeans that the DFT can be reduced to N multiplications of the data withthe eigenvalues. But as emphasized above, this observation is onlypartly true, i.e., although the performance gains from the diagonalrepresentation are indeed realized, these are offset by thecomputational cost of setting up the data in a transformed basis asillustrated above. But even with this in mind, certain applicationsspend most of the time in the forward and inverse transform part of thecalculations (see FIG. 2), so the main point is to explore applicationswhere the initial and final data transformations are needed only once,and then subsequent DFT calculations can proceed as simply Nmultiplications of the entries along the main diagonal. The performancegain for a single DFT calculation can then be approximately log(N) timesfaster. Assuming for simplicity a power of 2 for the data size(N=2^(n)), we can estimate a performance improvement for large n to be nlog(2)˜10¹, or about an order of magnitude faster than the equivalentFFT calculation, and for arbitrary frequency spacing.

FIG. 3 shows an example embodiment of the present invention and is thusdirected to an image-based wavefront sensing method of and means forfocusing by determining phase difference between received images,removing that phase difference by calculating a DFT for each signal anditeratively multiplying the resultant until a convergence is achieved byaccepting point-source stellar objects to recover embedded optical phaseinformation, receiving data from a sensor array wherein the dataincludes embedded phase information, calculating a similarity transformand pre-transforming the data once at the beginning of the calculation.

FIG. 3 further shows, in an example embodiment of the present invention,an image-based wavefront sensing method of and means for iterativelyprocessing the data using arbitrary frequency spacing in a diagonalfashion until a wavefront convergence is achieved to determine a phasedifference by implementing the DFT in a 1 dimensional linear array usingjust N operations, where 0<N<N×log(N)<∞, wherein N×log(N) is the numberof operations of an FFT counterpart, back-transforming the data at theend of the calculation in a single step and focusing the received imagesbased on the determined phase.

While the invention is described with reference to an exemplaryembodiment, it will be understood by those skilled in the art thatvarious changes may be made and equivalence may be substituted forelements thereof without departing from the scope of the invention. Inaddition, many modifications may be made to the teachings of theinvention to adapt to a particular situation without departing from thescope thereof. Therefore, it is intended that the invention not belimited the embodiments disclosed for carrying out this invention, butthat the invention includes all embodiments falling with the scope ofthe appended claims. Moreover, the use of the terms first, second, etc.does not denote any order of importance, but rather the terms first,second, etc. are used to distinguish one element from another.

1. An image-based wavefront sensing method of focusing by determiningphase difference between received images, removing that phase differenceby calculating a discrete Fourier transform (DFT) for each signal anditeratively multiplying the resultant until a convergence is achievedcomprising the steps of: accepting point-source images to recoverembedded phase information: receiving data from an external sourcewherein the data includes embedded phase information; calculating asimilarity transform; iteratively processing the data in a diagonalfashion until a wavefront convergence is achieved to determine a phasedifference; and focusing the received images based on the determinedphase.
 2. The wavefront sensing method of claim 1 wherein calculating asimilarity transform further includes pre-transforming the data once atthe beginning of the calculation.
 3. The wavefront sensing method ofclaim 2 wherein iteratively processing further includes a single step ofback-transforming the data at the end of the calculation,
 4. Thewavefront sensing method of claim 3 wherein said point source imagesinclude stellar objects.
 5. The wavefront sensing method of claim 4wherein said embedded phase information is of an optical nature.
 6. Thewavefront sensing method of claim 5 wherein said external sourceincludes a sensor array.
 7. The wavefront sensing method of claim 6wherein the step of iteratively processing further includes the step ofdetermining the phase difference by utilizing a discrete Fouriertransform (DFT).
 8. The wavefront sensing method of claim 7 includingimplementing said DFT in a 1 dimensional linear array using just Noperations, where 0<N<N×log(N)<∞, wherein N×log(N) is the number ofoperations of an FFT counterpart and N is an integer.
 9. The wavefrontsensing method of claim 8 including iteratively processing usingarbitrary frequency spacing.
 10. An image-based wavefront apparatus forfocusing by determining phase difference between received images,removing that phase difference by calculating a discrete Fouriertransform (DFT) for each signal and iteratively multiplying theresultant until a convergence is achieved comprising: means foraccepting point-source images to recover embedded phase information:means for receiving data from an external source wherein the dataincludes embedded phase information; means for calculating a similaritytransform; means for iteratively processing the data in a diagonalfashion until a wavefront convergence is achieved to determine a phasedifference; and means for focusing the received images based on thedetermined phase.
 11. The wavefront sensing apparatus of claim 10wherein said means for calculating a similarity transform furtherincludes means for pre-transforming the data once at the beginning ofthe calculation.
 12. The wavefront sensing apparatus of claim 11 whereinsaid means for iteratively processing further includes means for ofback-transforming the data at the end of the calculation in a singlestep,
 13. The wavefront sensing apparatus of claim 12 wherein said pointsource images include stellar objects.
 14. The wavefront sensingapparatus of claim 13 wherein said embedded phase information is of anoptical nature.
 15. The wavefront sensing apparatus of claim 14 whereinsaid external source includes a sensor array.
 16. The wavefront sensingapparatus of claim 15 wherein the means for iteratively processingfurther includes means for determining the phase difference by utilizinga discrete Fourier transform (DFT).
 17. The wavefront sensing apparatusof claim 16 further including means for implementing said DFT in a 1dimensional linear array using just N operations, where 0<N<N×log(N)<∞,wherein N×log(N) is the number of operations of an FFT counterpart. 18.The wavefront sensing apparatus of claim 17 further including means foriteratively processing using arbitrary frequency spacing.